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ABSTRACT 

This paper describes the development and testing of a general relativistic magneto- 
hydrodynamic (GRMHD) code to study ideal MHD in the fixed background of a Kerr 
black hole. The code is a direct extension of the hydrodynamic code of Hawley, Smarr, 
and Wilson, and uses Evans and Hawley constrained transport (CT) to evolve the mag- 
netic fields. Two categories of test cases were undertaken. A one dimensional version 
of the code (Minkowski metric) was used to verify code performance in the special rel- 
ativistic limit. The tests include Alfven wave propagation, fast and slow magnetosonic 
shocks, rarefaction waves, and both relativistic and non-relativistic shock tubes. A 
series of one- and two-dimensional tests were also carried out in the Kerr metric: mag- 
netized Bondi inflow, a magnetized inflow test due to Gammie, and two-dimensional 
magnetized constant-/ tori that are subject to the magnetorotational instability. 

Subject headings: Black holes - magnetohydrodynamics - instabilities - stars: accretion 



1. Introduction 



Accretion into black holes is believed to account for a wide variety of astrophysical phenomena, 
from solar-mass black holes in X-ray binaries, to supermassive black holes in active galactic nuclei. 
The basic theory of black hole accretion was laid out more than thirty years ago in a number of now- 
classic papers (e.g. Novikov & Thorne 1972; Shakura & Sunyaev 1973; Lynden-Bell and Pringle 
1974). Since then, theoretical progress has been steady, but perhaps at a slower rate than might 
have been hoped. The problem is complex, involving time-dependent, three-dimensional dynamics 
of magnetized plasmas in the relativistic potential of Kerr black holes. Solutions are impossible to 
obtain analytically beyond simplified cases that rely upon time-stationarity and spatial symmetry. 
Because of this, numerical experiments must play an increasingly important role in driving theory. 
The importance of numerical simulations in the investigation of such complex physics has long been 
recognized, but until recently the necessary computational hardware was not available. 
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We are now in a position to pursue full 3D simulations with general relativistic MHD codes. 
Again, the basics were laid out thirty years ago by Wilson (1972) who pioneered both hydrodynamic 
and magnctohydro dynamic (MHD) two-dimensional simulations of black hole accretion (Wilson 
1975; 1977; 1978). Since then there have been a significant number of efforts by a number of groups 
to develop accretion simulations for both Newtonian and relativistic gravitational fields, and for 
both hydrodynamics and magnetohydrodynamics. This paper details our current effort to develop 
a fully general relativistic MHD accretion code. 

Why focus on general relativistic MHD? In fact, for many applications Newtonian gravity 
(or a pseudo-Newtonian mock-up of a black hole) is sufficient. In some important cases, however, 
it is not. The detailed physics of gas interacting with a black hole provides, in principle, a rare 
test of strong-field general relativity, but only if we understand the unique properties of accretion 
in the fully relativistic case. Distinctive and astrophysically interesting effects are expected from 
the Kerr metric, including those due to the lack of a stellar surface, the presence of an innermost 
stable circular orbit, and the dragging of inertial frames. The need to include MHD in an accretion 
simulation is also now quite clear. Magnetic fields play an essential role in the outward transport 
of angular momentum in accretion disks through the action of the magnetorotational instability 
(MRI; Balbus & Hawley 1998). This, in turn, means the simulations must be three-dimensional, for 
only in three dimensions can there be a self-sustaining (and accretion-sustaining) MHD dynamo. 

In this paper we describe an algorithmic approach to the three-dimensional equations of general 
relativistic MHD that is based upon techniques first developed for axisymmetric hydrodynamics 
around black holes (Hawley, Smarr k, Wilson 1984a, 1984b; hereafter HSWa and HSWb). The hy- 
drodynamic portion of the code was recently tested and employed in three-dimensional simulations 
in the Kerr metric by De Villiers k, Hawley (2002). Here we cast the equations of relativistic MHD 
in a form that is compatible with the HSW formulation. We will describe how those equations are 
evolved numerically. We will also discuss a suite of tests suitable for relativistic codes, and show 
how our code performs against those tests. 

It is important to recognize that the development of numerical algorithms in any area of physics 
with few analytic solutions, and even fewer experimental observations, is likely to be difficult. 
Progress will be made only by exploring a variety of approaches. Such a situation necessarily 
precludes the development of only one numerical algorithm, or even a single code implementation 
of a promising algorithm. Some approaches will have strengths for certain applications, but it is 
unlikely that one numerical technique will be optimal for all problems. In any case, experience 
gained, both with what works and what does not, will be essential as we work toward developing 
a consensus within the community. With this in mind, our own efforts have concentrated on 
extending the established general relativistic hydrodynamic approach of HSW, since our main area 
of interest, that of accretion flows, has been well handled (in the hydrodynamic case) by this solver, 
and by similar algorithms designed to perform Newtonian MHD simulations. Our test suite has 
been developed in conjunction with Gammie, McKinney, &; Toth (2002), who use an alternative 
scheme for general relativistic MHD. 
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The plan of this paper is as follows. In §2 we write the equations of general relativistic MHD 

in the form that they will be solved in the code. In §3 we describe the numerical implementation of 
our algorithm. In §4 we describe the one-dimensional Minkowski (special relativistic) test problems. 
In §5 we describe the test problems in the Kerr and Schwarzschild metrics, and present our test 
results. The summary and conclusions are presented in §6. 



2. Equations of General Relativistic Magneto-Hydrodynamics 

We wish to study the evolution of a magnetized fluid in the background spacetime metric of a 
Kerr black hole. We work in Boyer-Lindquist coordinates, for which the line element has the form, 

ds^ = gtt dt^ + 2 gt(i) dt dcf) + grr dr^ + gee dO"^ + g^^ dcfP'. (1) 

In keeping with Misner, Thorne, & Wheeler (1973), we use the metric signature (—,+,+,+), 
along with geometrodynamic units where G = c = 1; the black hole mass is unity, M = 1. The 
determinant of the 4-metric is g, and \^—g = «\/7 where a is the lapse function, a = I/a/— 5**, 
and 7 is the determinant of the spatial 3-metric. 

The state of the relativistic test fluid at each point in the spacetime is described by its density, 
p, internal energy, e, 4-velocity U^, and isotropic pressure, P, which is related to the first two 
scalars via the equation of state of an ideal gas, P = pe (7 — 1), where 7 is the adiabatic exponent. 

The equations of general relativistic MHD are the law of baryon conservation, 

V^(pC/^) = 0, (2) 

where is the covariant derivative, the conservation of stress-energy, 

V^.T^"' = 0, (3) 

where Tf^" is the energy-momentum tensor for the fluid, and Maxwell's equations, 

Vf.Ff"' = A7rr, (4) 
V/F'^'^ = 0, (5) 

where F^" is the electromagnetic field strength tensor, Ji^ = {pc, J*) is the current 4- vector, and pc 
the charge density. The dual tensor is defined as 

^i^'^" = ^e'^'^^^Fs^, (6) 

where e*^"^'^ = —(1/-^/^) [p,uSj] is the contra-variant form of the Levi-Civita tensor. Maxwell's 
equations are supplemented by the equation of charge conservation J'* = 0. 
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The energy- momentum tensor consists of perfect fluid and electromagnetic parts, T^'^ = 
T^/iuid) + Tj^EMY and has the form 



1 



Ti''' = phU^'U'' + Pg^''' + — [Ft" F'''' - -F^^ F'^P g^""" j . (7) 

Following Lichnerowicz (1967) and Anile (1989), we define the magnetic induction and the electric 
field in the rest frame of the fluid, 

B^' = *F''''U^, (8) 
= F'^'U^. (9) 

Wc adopt the ideal MHD limit and assume infinite conductivity (the flux- freezing condition), 
wherein the electric field in the fluid rest frame is zero, i.e., Fi^'^ = 0. 

We combine the definition of the magnetic induction (8) with (6) and the condition for infinite 
conductivity to obtain 

i^/.. = e„/3;..5"i7^, (10) 

where e^,/^^ 

= \r-9\M'^^l\- The orthogonality condition 

B^'U^ = Q (11) 

follows directly from (8) and the anti-symmetry of F^^ . 

Using these results, it is possible to rewrite the electro-magnetic portion of the energy-momentum 
tensor as 

Ti;^M) = {^1^9'''M? + U^U^m?-h^h^y (12) 

where = BI'/Va^, and \\bf = g'^'b^b^. 

The induction equation (5) is solved using the equivalent form 

dsF^p + daFf^s + dpFsa = 0. (13) 

Working directly with F^j^ replaces covariant derivatives with simple coordinate derivatives. This 
is the basis for the Constrained Transport (CT) formulation of Evans &; Hawley (1988; hereafter 
EH88) for the induction equation. In CT we make the identification 

= F^e , = Fr^ , B'^ = Fer, (14) 

where arc the CT magnetic field variables. With t for one of the indices in (13) we get the evo- 
lution equations for B^ , while permuting over the spatial indices yields dr F^q + Oq Fj.^ + Fg^. = 0, 
which is the familiar divergence-free condition. The CT magnetic field is considered the fundamen- 
tal expression of the magnetic field for evolution in the induction equation. Equation (10) relates 
the CT field to the field B^^ (or b^^) which is fundamental in the definition of T^^M)- 
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The induction equation (5) can also be rewritten by substituting definitions, 

Va - 6" t/^) = 0. 



(15) 



By expanding this equation using the product rule and applying the orthogonality condition (11), 
we obtain the identity 

U^h^'V^U'' = Q. (16) 

Proofs of the results for (10), (12), and (15) can be found in Appendix A. 

In preparation for discretizing these equations we list the specific variables that will be used 
in the code. First, the transport velocity (also known as the coordinate velocity) is defined 



= — 

where = W/a, and W is the relativistic gamma-factor, which can be expressed as 



W 



5" [att + 29t^ + grr {y'f + gee iy' f + 9h (^'^)') 



-1/2 



(17) 



(18) 



a result which follows directly from the definition of transport velocity and the normalization 
condition on the 4-velocity, 

g^,U>'U'' = -\. (19) 

The magnetic field of the fluid is described by two sets of variables, the constrained transport 
(CT) magnetic field, ^ and magnetic field 4-vector . The latter is fundamental to the definition 
of the total four momentum, 

S^, = {ph -r\\hf)WU^, (20) 

and the normalization condition 

1 2 



g^''S^S, = - {ph + 



(21) 



which is algebraically equivalent to (19). Finally, we define auxiliary density and energy functions 
D = pW and E = D e. The set of variables D, E, 5^, B*, V^, and bf^ will be those on which the 
numerical scheme is built. 

Since we are working in a coordinate basis, we also record here the identities for the divergence 
of a four- vector and a tensor 



1 



1 



(22) 
(23) 
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2.1. The Equations of Baryon, Energy, and Momentum Conservation 

The equation of baryon conservation (2) is unaltered by the presence of magnetic fields, and can 
be expanded in terms of the code variables to read (keeping in mind that the metric is stationary) , 

dtD + ^ dj {D y^) = 0. (24) 

The equation of energy conservation is derived by contracting (3) with Ui,, 

f/,V/,r'^" = C/,V/,|(p/i+||6||') f/^f/-+ ||p+i^^5'^--6'^6-| =0. (25) 

By using identity (16) and the law of baryon conservation (2), we recover the local energy conser- 
vation law 

V^,{peU^') + PV^U^' = {), (26) 

which is unaffected by the presence of magnetic fields in the ideal MHD limit and corresponds to 
an equation of entropy conservation. Applying the definition for the auxiliary energy function E, 
the energy equation is rewritten as follows: 

dt{E) + ^di{^EV')+Pdt{W) + ^di{^WV')=Q, (27) 

which is the same as in HSW. 

The momentum conservation equations follow from 

T\ = I (p + \\hf) U,+ (^P+ ,5^ - ft'* 6, 1 = 0. (28) 

Using the definition of momentum Si, the first term in the preceding expression can be rewritten as 
SuV^^/a and simplified to SyS^^/aS*, and also using T^j^TI^^ = ^IP g^y = —^U^nd^g^"', 
which holds for any symmetric tensor, we rewrite the momentum equation as 

4= a. V7 s„ v" + ^ ^ a. /"^ + a„ (p + H! J _ _i_ a„ „ ^ 1„ - i 6„ 6, a„ = 0. 

(29) 

To obtain the final form of the equations, multiply (29) by the lapse a, split the index into its 
space {i) and time {t) components, and restrict u to the spatial indices (j) only: 

dt{Sj-abjb') + ^diVl {SjV'-abjb') + ^ (^^-ab^b,^ dj g^^ + adj (^P + = 0. 

(30) 

The u index can be restricted to the spatial indices because the equation that arises from u = t 
for the time components of momentum and magnetic fields is redundant, corresponding to a total 
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energy conservation equation. In our formalism, we solve equation (27) separately for the internal 
energy. The time component of is obtained from the normalization condition on the momentum 
(21), and on the relations between components of the magnetic field. 

We note that there are alternate ways to write this equation that, while analytically equivalent, 
correspond to significant differences in numerical implementation. One is to evolve directly the total 
stress energy T^'^ and solve algebraically for the primitive variables. This approach is taken by 
Koide, Shibata, k Kudoh (1999), Komissarov (1999), and Gammie, McKinney & Toth (2002). Yet 
another alternative is to write the divergence of the electromagnetic stress-energy tensor in the 
form JfxF^^^ (the familiar J x B force from Newtonian MHD) and regard it as a source term for 
the hydrodynamic momentum evolution; this approach was used by Wilson (1978). 



2.2. The Induction Equation 

The equations of energy and momentum conservation have been expressed in terms of the 
magnetic field 4-vector 6'^ since this yields expressions structurally compatible with the existing 
hydrodynamic equations. This has important advantages in the numerical implementation of these 
equations. However, the induction equation is evolved with the constrained transport (CT) ap- 
proach of EH88 using the variables B^. The practical consequence is that we must translate between 
the two sets of variables. 

Working with the second of Maxwell's equations, (5), we have previously identified the space- 
space components of F^i, with the CT field B^. It is easy to show that we obtain the same result 
using (15), 

(U" - b" U^) = W y - b" V^") = 0. (31) 
By splitting the above equation into two pieces and defining the CT magnetic field as 

B' = VA^^W{b' -V'b^), (32) 

it is possible to write 

dj{B^) =0 {iy = 0), (33) 

dt [B') - dj {V' B^ -B'V^) =0 {iy = i), (34) 

the starting point for the EH88 CT scheme. In addition, the time-space components of F^i/ can be 
obtained by applying the condition of infinite conductivity, F>^'^ Uy = 0, to obtain 

Ftr = V^B^ - V^B^,Fte = V" B'^ - B\ Ft4, = B' - VB'^. (35) 

These can be identified with £^ , the CT electromotive forces (EMFs, EH88). 

The divergence-free character of the magnetic field is maintained by evolving the induction 
equation using these CT variables. Hence the magnetic induction b^ is derived from the CT 
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magnetic field. We obtain the spatial components using (32). To obtain 6*, we use the definition 
of the magnetic induction (8), 

b'' = ^*F'^-C/, = ^—ef^-'^U^Fs^, (36) 
V47r 2 V 47r 

and we expand this definition and substitute for the Faraday tensor. This allows us to obtain the 
time-component of the magnetic field, 



V47ra2^ V ^'"'^ gtt g<t><l> - {gt<t>y 

This relationship can also be obtained from the orthogonality condition (11). Using this result and 
the relations in (32) we obtain, after some algebra, a useful expression for 

11,^11 ^ _i_fg:!^g!^ fg^' \ (38) 

A-K-fW^ ygrr gee g^^ g'^"*' - {g^'f') ) 

1 fvB'' B^ (c/** -g^^)B^V 
^47r7a2 y'-'- ^ gO(> ^ g** g^^ - {g^^ f J ' 

The advantage of this form is that it guarantees that \\b'^\\ is positive. The derivation is supplied 
in Appendix A. 



3. Implementing GRMHD 

The GRMHD code evolves time-explicit, operator-split, finite difference forms of equations 
(24), (27), (30), and (34). Variables are placed on a fixed spatial grid using Boyer-Lindquist 
coordinates. They are advanced in time over a timestep size At that remains fixed for the duration 
of the simulation, and is determined by the extremal light-crossing time for a grid zone, as described 
in HSWb. The evolution algorithm is a three-dimensional generalization of the solver described in 
HSWb extended to include the contribution of magnetic fields. The hydrodynamic portion of the 
solver was described in De Villiers and Hawley (2002). 

One timestep consists of three sequential sub-steps: the induction step, the transport step, 
and finally the source step. The induction step updates the magnetic field vector B^ using the 
Constrained Transport framework of EH88. This is discussed in greater detail below. After the 
update, the CT field is transformed to the magnetic field in the fluid rest frame, 6*. Velocities 
are computed using the normalization condition during the source step. 

Three versions of GRMHD have been used in development and testing: (1) a ID version with 
Minkowski metric to do Alfven and shock wave tests, described here in §4, (2) a 2D axisymmetric 
version in the Kerr metric to perform the tests described in §5, and (3) the full 3D production 
version of the code. The 3D version of the code uses message passing parallelism with domain 
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decomposition, where the global grid is partitioned into subgrids, with each subgrid assigned to a 
processor. Data on each subgrid are evolved independently during a timestep and data on subgrid 
boundaries are exchanged when required through message-passing calls. This results in a highly 
scalable code that exhibits good speedup over the full range of practically realizable subgrids. 



3.1. Transport and Source Steps 



The transport step, as its name implies, handles the transport terms in the continuity, energy, 
and momentum equations, 



[dtX] 



1 



transport 



V7 



(39) 



Other than the modified definition of the inertia, these expressions are completely analogous to 
those in HSWb. 



The source step handles the remaining terms in the energy and momentum equations. The 
energy terms 

[dtE]source = -P dtW - [VlWY^ (40) 



,^ = -PdtW -—di y^wv'] 

are evolved as described in HSWb. The momentum source terms contain new elements, namely 



[dtSjlource = -<^dt[bjb'] + 



V7 



2 ^ 



P + 



■ (41) 



In the following expressions, superscript indices for the current time step are understood. The 
spatial differencing of some value / is shown in compact form through the shift operator Di(f); any 
averaging of zone- and face-centered quantities that may be required to properly center calculations 
is done as described in HSWb. 

The time derivative of the magnetic fields is evaluated from the newly-obtained values from 
the induction step, and stored values from the previous time step. 



Sr' = Sj + a[{b,b'Y''^'^-{b,b') 
The magnetic gradient term discretizes readily, 

At Di{^abjb') 



= Sj + 



V7 



(42) 



(43) 



as does the pressure gradient term, 



= Sj-Ata 



Dj[p+\\bf/2) 



Ax^ 



(44) 
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Finally, the accelerations due to gradients in the metric, 




At \SeS, 



-abf^be dj g' 



(45) 



2 St 



are computed using 5* and St obtained from Sj using momentum normaHzation (21); the metric 
derivatives, djg^^^, are evaluated analytically over the whole grid and stored in memory at the 
beginning of the simulation. 

In the source step, after the pressure gradient and acceleration terms have been computed, 
the transport velocity = S'^/S^, and the Lorentz factor W = aU* are obtained using velocity 
normalization. To guarantee that the normalization is well-behaved in evacuated zones, we work 
with the momenta Si rather than applying the normalization condition directly to the velocities. 
Although the two normalizations are analytically equivalent, use of the momentum protects against 
numerical rounding errors that cause spurious floating point exceptions when extremely tenuous 
material is moving at velocities close to that of light. 

The source step also includes an artificial viscosity in both the energy and momentum equations 
to provide a mechanism to increase the entropy of the fluid in shocks. The artificial viscosity is 
unchanged from HSWb, except for modifying the definition of the inertia to include the magnetic 
energy. 

3.2. The Induction Step: A General Relativistic Form of the MOCCT Algorithm 

The CT framework of EH88 is based on two basic ideas: the use of F^y as the fundamental 
variable so that the induction equation reduces to simple partial derivatives, and a staggered-grid 
centering of terms so that the numerical divergence of the CT magnetic field is constrained to be 
zero. The CT magnetic field component is located on the j face of a zone (a cube in three 
dimensions), and the CT EMFs, £'^, that are differenced to evolve are located on the zone edges 
that define the j face. 

The CT formalism imposes no stringent conditions on exactly how the EMFs are computed. 
EH88 noted that some of the magnetic field variables appear in the discretrized equations for the 
EMFs in the guise of transport terms, requiring that they be upwinded for stability, while others 

appeared as shear terms, apparently requiring no special treatment. EH88 used simple space- 
centered differencing for those terms, but Stone &: Norman (1992) found that this prescription is 
inadequate. The results from simple Alfven pulse propagation tests where the Alfvenic velocity 
exceeds that of the background fluid found substantial undamped dispersion error, resulting in 
excessive zone-to-zone noise. 

Stone & Norman (1992) proposed to calculate the CT EMFs by solving a one dimensional 
analytic linear Alfven wave equation to obtain improved values of B^ and in the construction of 
the EMFs. This approach became known as the Method of Characteristics Constrained Transport 
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(MOCCT). Since then several variants have been proposed (e.g. Hawley & Stone 1995; Clarke 

1996) to further improve the computation of Newtonian MHD. While the traditional approach to 
MoC does not carry over readily to GR MHD, it is nevertheless possible to address the problem 
by using the same motivation that underlies MOCCT. Instead of using an analytic Alfven wave 
solution we use a numerical solution of the simplified one-dimensional induction and momentum 
equations to obtain the improved EMFs. 

Our implementation of CT for the GRMHD code is accomplished as follows. Following the 
example of Hawley & Stone (1995), the starting point for constructing, for instance, the z-EMF, 

= - B'' , (46) 

lies in obtaining the discretized solutions to the ID incompressible MHD equations. On a staggered 
mesh, is located on zone edges in such a way that it is constructed from x- and y-components 
relocated in the y- and .x-dircctions respectively. To take a specific example, we would like obtain 
an edge-centered estimate B^* to construct the above EMF. We do this by solving the ID linearized 
induction equation in the x-direction (i.e. dy = dz ^ 0) for small transverse (y) perturbations. This 
simplified induction equation has the form, 

= 0, (47) 
dt By = 5^ {yy) - 9^ [B^) . (48) 

By applying the usual differencing we obtain for the left-hand side, 

^tB' = ^^, (49) 

where the overbar denotes the original EH88 upwinding. We can now readily obtain estimates of 
the starred quantity, 

By*=m+l^ [-v^ D^isy) + B^D^(yy)] , (50) 

2 Ax 

where the factor of 1/2 ensures proper time-centering of the starred quantity (the starred values 
are centered between the n and n + 1 time levels). A similar approach can now be used to generate 
B^*, and produce an improved estimate of S^. The procedure is just as easily applied to the other 
two EMFs to complete the process for all CT variables. Note that because of the fully covariant 
nature of the CT equations, this process applies to any coordinate system. 

For consistency, we also need to evolve the linearized momentum equation for the transverse 
velocities. Care must be taken in deriving these equations since, in contrast to the CT induction 
equation, metric terms are explicitly involved in the expressions. Continuing with the example for 
and regarding the x direction to be either r or 9, and y a direction transverse to this, we obtain 



{dt + c d^) yy-K B'^d^ By o. 



(51) 
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where, 



e = 



phW^ - 
phW^ + 



|2 ' 



(52) 
(53) 



(The details are found in Appendix A.) Note that the above equations are similar to the Newtonian 
results, but with the relativistic factors ^ and k. By simply applying the usual differencing, however, 
we readily obtain expressions for an improved estimate of the transverse velocity. 



dfVy 



yy* _ yy 
At 



where the overbar denotes the edge-centered average, from which, 

yy* = w + i ^ [-^ D,{vy) + kB'^d, By] . 



(54) 



(55) 



In the Kerr metric, the calculation in the (f) direction is somewhat different due to the g^'^ 
metric term. In this case we obtain 



5t+ (CF^ + A) a, 
where ^ and n are unchanged and 



yy - K B^dA,By ^ o. 



A 



2||6f fl*-^ 



gtt(^phW^ + \\bfy 



(56) 



(57) 



Combined with the starred quantities for the CT variables, we now have a complete set of 
variables from which to construct the EMFs. 



4. ID Test Cases - Minkowski Spacetime 

The first series of tests is designed to verify code performance in the special relativistic limit. 
For this purpose we use a one-dimensional version of the code and the Minkowski metric. The tests 
include Alfven wave propagation, fast and slow magnetosonic shocks, rarefaction waves, and both 
relativistic and non-relativistic shock tubes. Many of these tests have appeared previously in the 
literature in one form or another. In particular, Komissarov (1999) (hereafter referred to as K99) 
presented a series of challenging test problems for a special relativistic MHD Godunov scheme, and 
Stone et al. (1992) presented a test suite for Newtonian MHD codes. Our presentation will conform 
to the formats of these published results to facilitate comparison. 
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4.1. Alfven Wave Propagation 

The propagation of a linear Alfven wave is simple in conception, but, as demonstrated by 
Stone Sz Norman (1992), can be quite revealing as a numerical test. The special relativistic Alfven 
wave includes the displacement current normally neglected in Newtonian MHD. This is what limits 
wave speeds to less than c. 

We consider a fixed background magnetic field in Minkowski spacetime with constant fluid 
velocity and a small transverse perturbation with velocity and CT magnetic field B^. The 

induction and momentum equations for the perturbed quantities have the form given in equations 
(48) and (51). Equation (51) can be combined with the induction equation (48) by taking partial 
derivatives [{dt + dx) for the former and dx for the latter] to get an equation for , 

a| + (1 + dt dx + {V^f - n {B^f) dl\ yy = O (58) 
which describes the propagation of Alfven waves with speeds 



1 + ^^ ~\ ph + 



VA^^^ = ^-^:LV±_L^ ^ Mi^^^AforV^ = 0), (59) 



where r/^ = ||6||^/(/9/i H^^). (The equation for B"^ can be obtained in a similar manner.) Figure 
1 plots VA^^^ as a function of the background velocity of the fluid (V^) for four values of (3 and 
illustrates several key properties of the wave speeds: (1) for = 0, waves move with equal and 
opposite speed (for all (3); (2) wa^^^ is constrained to the upper bound of c = 1; (3) for a given /?, 
there is a Vq that yields a stationary wave, ■ua'^"-' = 0; (4) for > Vq, both waves move in the 
same direction, with ■ua^"-' lagging behind the main fluid flow, and va^'^^ leading; (5) as —>■ c, 
all wave speeds converge, va^"^"^ — *■ c. 




Fig. 1.— Plot of va'^^^ for /3 = 10'^ (solid line), /3 = 10"^ (short dashes), /? = 10"^ (dash-dot), 
P = 1 (long dashes). 

We describe four tests, listed in Table 1, carried out in a periodic box of width Xmax — Xjjiin — 
3.0. In each case the fluid has constant density p = 1 and energy density e = 10~^. Individual tests 
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are distinguished by the amphtude of the x-component of the CT magnetic field B^, given here in 
terms oi (3 = ^J2P/\\b\\^ , and the amphtude of the background x-component fluid velocity . The 
waves are initialized with a square pulse in the transverse velocity with amplitude Aq = 10~^, 
namely 

r ^0 for 1.0 < X < 1.5 

\ -Aq for 1.5 < X < 2.0 ■ ^ ' 

The transverse magnetic field is initially zero, = 0. 

Table 1 summarizes the results and Figure 2 overlays for each test at time t final at three 
grid resolutions, 512, 1024, and 2048 zones. In test ALF2 the fast and slow pulses have wrapped 
around and re-entered the grid from the left, the fast pulse having done so twice. In ALF4, the 
fast pulse has also wrapped around while the slow pulse is stationary. (In this linear regime the 
pulses overlap cleanly provided that the spatial resolution is sufficiently high to minimize losses 
due to numerical diffusion.) Measuring from the center of each pulse (originally at x = 1.5) we 
verify that the observed pulse speeds agree with (59). With a stationary background, ALFl, the 
pulses separate from the initial perturbation with equal and opposite velocities and also with equal 
amplitudes such that = + In ^ moving background, the splitting of the velocity pulses 
is no longer symmetrical, as is especially obvious in ALF2 (although = A^ + A~ still holds). 
The slow velocity pulse has a greater amplitude than the fast pulse. We find that the ratio of the 
velocity amplitudes scales inversely with their Lorentz factors, 

^ = (61) 

A- W{vA+) ^ ' 

It is easy to verify that this holds for each model in Table 1. The pulses in the magnetic field 
variable, , on the other hand, are equal in magnitude and opposite in phase for all V'^ . The 
numerical values given in Table 1 can be verified by direct computation using results presented in 
Appendix B. 

Table 1: Results of ID Alfven Pulse Tests. 



Model 


yx 






(xio-«) 


w+ 


va" 


A (xio-*) 


w- 


i final 


\\By\\ (xio-3) 


ALFl 


0.000 


0.001 


+0.96 


5.00 


3.76 


-0.96 


5.00 


3.76 


0.9 


6.714 


ALF2 


0.800 


0.100 


+0.90 


3.64 


2.26 


+0.63 


6.36 


1.29 


5.8 


4.888 


ALF3 


0.100 


0.010 


+0.79 


4.62 


1.64 


-0.71 


5.38 


1.41 


1.1 


2.729 


ALF4 


0.200 


0.315 


+0.38 


4.80 


1.08 


0.00 


5.20 


1.00 


6.2 


1.897 



4.2. Magnetosonic Shock and Rarefaction Tests 



Shock waves and rarefactions test a code's ability to respond to discontinuities and to maintain 
the jump conditions that correspond to the conservation properties of the flow. In this section we 
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Fig. 2. — Transverse velocity at time t final for Alfven wave tests (from left to right) ALFl, 
ALF2, ALF3, and ALF4. Three different resolutions are overlaid: 512, 1024, and 2048 zones. 



consider a number of one-dimensional shock tests, both relativistic and non-relativistic, that involve 
nonzero transverse magnetic fields. 

We begin with relativistic magnetosonic shock wave tests; the initial left {x < 0.0) and right 
{x > 0.0) states are given in Table 2. In each case we verified that these states satisfy the relativistic 
shock invariants (see Appendix C). In this and in subsequent tables we note those tests that were 
described in K99.^ 

The first test, the slow magnetosonic shock, is shown in Figure 3, which plots U^, p, gas 
pressure P, and magnetic pressure Pmag at t = 2.0, by which point the shock front has moved from 
X = 0.0 to a; = 0.5. These graphs are quite similar to the results of K99. There is a small step in 
density in the high-density medium to the right of the moving shock. When we repeat our tests at 
lower resolution, our step feature resembles the "wiggle" visible in K99. The origin of this small 
step may be due to a small discrepancies in the initial state (i.e. the data in K99, which we use 
here, differ slightly from what we compute using the method outlined in Appendix C). 



2.0 - 
1.5 ^ 
1.0 - 
0.5 : 
0.0 _ 



Fig. 3. — Plot of (left), Pgas (bottom curve, center) and Pmag (top curve, center), and p (right) 
at t = 2.00 for the K99 Slow Magnetosonic Shock. 

Figure 4 shows plots of U^, p, gas pressure P, and magnetic pressure Pmag for the fast magne- 
tosonic shocks at t = 2.5. In the case of the stationary shock (I), we see that artificial viscosity has 
slightly altered the shock front, and that there is an overvalue in the left state [/^, and undervalue 
in p and Pgas- In the case of moving shock (II), the shock front has traveled to x = 0.5, as expected 
for the given shock speed and the value of t final, and the left and right states agree with the analytic 



^K99 magnetic field variables B}fgg are related to our CT- variables by a normalization factor, = v'JttjBJ, 
Also, a revised version of K99 Table 2 was given by Komissarov (2002). 
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Table 2: Initial States for ID Magnetosonic Shock Tests. 



Test 


Left State 


Right State 


Grid 


t final 


Slow Shock 




= (1.53,0.0,0.0) 


[/' 


= (0.9570, -0.6820,0.0) 


512 


2.0 


(V = 0.5) 




= (0.838, 0.0, 0.0) 


V 


= (0.620, -0.442, 0.0) 






(K99) 




= v47r(10.0, 18.28, 0.0) 


B 


= v47r(10.0, 14.49, 0.0) 








P-- 


= 10.0 p = 1.0 


P- 


= 55.33 p = 3.322 






Fast Shock I 


{/* 


= (1.780,0.114,0.0) 


W 


= (0.922, 0.403, 0.0) 


1024 


2.5 


(y = 0.0) 




= (0.870,0.056,0.0) 




= (0.650,0.284,0.0) 






{W « 2.0) 




= v4vr(3.33, 2.50, 0.0) 


JO? 

B 


= v47r(3.33, 4.52, 0.0) 








F - 


= 2.015 p = 1.406 


P - 


= 5.135 p = 2.714 






Fast Shock II 


U' 


= (1.780,0.114,0.0) 


W 


= (1.479,0.280,0.0) 


1024 


2.5 


{V = 0.2) 


V' 


= (0.870,0.056,0.0) 




= (0.818,0.155,0.0) 






{W ^ 2.0) 




= V47f(3.33, 2.50, 0.0) 




= \/4¥(3.33, 3.25, 0.0) 








P -- 


= 2.015 p = 1.406 


p -- 


= 2.655 p = 1.725 






Fast Shock III 


U' 


= (3.649,0.114,0.0) 


w 


= (0.715,0.231,0.0) 


128 


2.5 


{V = 0.2) 


V' 


= (0.964,0.030,0.0) 




= (0.542,0.185,0.0) 






{W ^ 3.8) 




= ^/4^(3.33, 2.50, 0.0) 




= \/4¥(3.33, 6.52, 0.0) 








p -- 


= 2.015 p = 1.406 


p -- 


= 34.99 p = 8.742 







values. In the case of the W ~ 3.8 shock (III), a small mismatch has developed in the post-shock 
gas pressure near the shock edge. In addition, the speed of the shock slightly exceeds the expected 
value 0.2. As discussed by HSWb and studied in detail by Norman &: Winkler (1986), the artificial 
viscosity used in the code to increase the entropy of the fluid through a shock produces sufficient 
heating for 1^ < 2, but increasingly underestimates the heating as W increases. The results of this 
shock test are consistent with these earlier findings. 

We next consider a pair of challenging rarefaction tests designed by K99. The initial left and 

right states of these tests arc given in Table 3. Figure 5 shows a plot of U^, gas pressure P and 
density p for the switch-off fast rarefaction at t = 1.0 and the switch-on slow rarefaction at t = 2.0. 
Our results generally agree with those of K99. The fast rarefaction extends from x = —0.5 to 
X = +0.6 and the slow rarefaction extends from x = —0.5 to x = -1-0.9. There is high-frequency 
noise in the plot of and P to the right of the slow rarefaction. We attribute this noise to 
an artifact of the routine that converts CT- variables to fluid-frame magnetic field. This test is 
particularly sensitive to the averaging techniques that are an unavoidable part of staggered-mesh 
discretizations. This noise is reduced at lower grid resolution, where increased numerical diffusion 
acts to damp out the oscillations. 




Fig. 4. — Plot of U^, Pgas and Pmag, and p at t = 2.50, for the Fast Magnetosonic Shock I (left 
column); the Fast Magnetosonic Shock II (center column); and the Fast Magnetosonic Shock III 
(right column). 



Tabic 3: Initial States for ID Rarefaction Tests. 



Test 



Left State 



Right State 



Grid tfinal 



Switch-off Fast 

Rarefaction 

(K99) 



U' = (-2.0,0.0,0.0) 
= (-0.894,0.0,0.0) 
= ^/4^(2.0, 0.0, 0.0) 

p = 1.0 p = 0.1 



P = 



-- (-0.212,-0.590,0.0) 
= (-.180,0.500,0.0) 
: \/4^(2.0, 4.71, 0.0) 
10.0 p = 0.562 



2048 1.0 



Switch-on Slow 

Rarefaction 
(K99) 



U' = (-0.765,-1.386,0.0) 

V' = (-0.409,-0.740,0.0) 
B' = ^47(1.0,1.022,0.0) 
p = 0.1 p= 1.78 X 10-3 



U' = (0.0,0.0,0.0) 

V' = (0.0,0.0,0.0) 

B' = \/47(1.0,0.0,0.0) 

p = 1.0 p = 0.01 



2048 2.0 



Next we consider a series of shock tube tests that combine strong shocks and rarefactions of 
various types. We ran two test problems from K99, the Brio & Wu (1988) shock tube test as 
described in Stone et al. (1992), an alternate version of the Brio &; Wu shock tube using the same 
initial state but with 7 = 4/3 instead of 7 = 2 to ensure that the sound speed is not relativistic, 
and a non-relativistic shock tube from Ryu & Jones (1995). In the non-relativistic shock tubes 
the values of pressure and density are rescaled to ensure that the sound speed is much less than c. 
Table 4 summarizes the various states for the shock tube tests, following the convention of Stone 
et al. (1992). Figures 6, 7, 8, and 9 show our numerical results for the evolution of the shock tubes 
at a grid resolution of 2048 zones. Some ringing is present in the non-relativistic shock tubes, at 
the location of the compound wave and at the slow shock in the density variable, for instance, 
whereas this ringing is absent from the relativistic shock tubes. This appears to be a manifestation 
of decreased dispersion error at relativitistic velocities. 
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Fig. 5. — Plot of U^, Pgas, and p at t = 1.00 in the K99 Fast Rarefaction Test (left column); and 
plot of ?7^, Pgas, and p at t = 2.00 in the K99 Slow Rarefaction Test (right column). 



Table 4: Initial and Intermediate States for Shock Tube Tests. (7 = 4/3 unless indicated.) 



Test 


Variable 




Left 


FR 


SC 




CDr 


FR 


Right 


Shock Tube 1 


P 




1.00 


0.08 




0.85 


0.10 




0.10 


(K99) 


P 




1000 


31.8 










1.00 


(yy = By = 0) 






0.00 


0.90 










0.00 


Shock Tube 2 


P 




1.00 


0.24 




0.63 


0.10 




0.10 


(K99) 


P 




30.0 


4.50 




16.5 


1.00 




1.00 


{yy = 0) 






0.00 


0.85 










0.00 








20.0 


9.14 










0.00 


Brio and Wu 


p 




1.00 


0.51 


0.68 


0.55 


0.36 


0.11 


0.13 


(Rel.) 


p 




1.00 


0.41 


0.60 


0.45 


0.45 


0.08 


0.10 




yX 




0.00 


0.40 


0.26 


0.28 


0.28 


-0.12 


0.00 




yV 




0.00 


-0.15 


-0.46 


-0.58 


-0.07 


-0.07 


0.00 




Byj^ 




1.00 


-0.14 


-0.37 


-0.42 


-0.42 


-0.83 


-1.00 


Brio and Wu 


p 




1.00 


0.67 


0.86 


0.71 


0.25 


0.12 


0.13 


(S92) 


P(xlO-4) 




1.00 


0.45 


0.73 


0.50 


0.50 


0.09 


0.10 


7 = 2.0 


t>^(xl0-3) 




0.00 


6.54 


4.51 


6.02 


6.02 


-2.76 


0.00 




t;2'(xlO-2) 




0.00 


-0.24 


-1.15 


-1.58 


-1.58 


-0.20 


0.00 




H2'(xlO-2)/\ 


/47r 


1.00 


0.57 


-0.37 


-0.54 


-0.54 


-0.89 


-1.00 


Ryu and Jones 


P 




1.00 


0.68 


0.74 


0.91 


0.61 


0.32 


0.40 


(RJ95) 


P(xlO-4) 




1.00 


0.53 


0.60 


0.86 


0.86 


0.28 


0.40 


7 = 5/3 


z;^(xl0-3) 




0.00 


6.66 


4.62 


3.03 


3.03 


-5.61 


0.00 




z;?'(xlO-2) 




0.00 


-0.90 


-1.48 


-1.24 


-1.24 


-0.80 


0.00 




^2^(xlO-2)/\ 


/47r 


1.00 


0.28 


-0.34 


-0.24 


-0.24 


-0.44 


-1.00 
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Fig. 6.— Plot of J7^, Pgas (and Pmag at right), and p aX t = 1.00 for K99 Shock Tube 1 (left 
column) and K99 Shock Tube 2 (right column). 




Fig. 7. — Plot of p (top row, left), Pgas (top row, center), /\/4lTT (top row, right), (bottom 
row, left), and (bottom row, center) at t = 1.00 for the Relativistic Brio and Wu Shock Tube. 
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Fig. 8. — Plot of variables p (top row, left), Pgas (top row, center), H^/\/47r (top row, right), 
(bottom row, left), and (bottom row, center) at f = 30.00, in the 7 = 2 Non-Relativistic Brio 
and Wu Shock Tube. 




Fig. 9. — Plot of variables p (top row, left), Pgas (top row, center), j^f^ (top row, right), 
(bottom row, left), and V"^ (bottom row, center) at t = 45.00, Non-Relativistic Ryu and Jones 
Shock Tube. 
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5. ID and 2D Test Cases - Kerr Spacetime 

To validate the general relativistic aspects of the code we require code tests that involve mag- 
netic fields and the Kerr geometry. Test problems with analytic solutions include one dimensional 
magnetized Bondi inflow, and a two-dimensional Bondi flow along a split monopole field; in these 
tests the field does not affect the dynamics of the problem. We also perform a onc-dimensional 
magnetized inflow test problem where the field does act dynamically (Gammie 1999). Finally, we 
have considered two-dimensional magnetized constant-/ tori that are subject to the magnetorota- 
tional instability. While there is no analytic solution to this problem, the results can be compared 
to those obtained in non-relativistic simulations. 



5.1. Magnetized Bondi Flow and Split Monopole Tests 

The analytic expressions for the hydrodynamic Bondi solution of HSWa are unchanged in 
the presence of a radial magnetic field (see proof in Appendix A). However, this result does not 
necessarily hold for a numerical solution to the Bondi problem. A non-trivial test of the magnetic 
components of the code consists, therefore, in maintaining the equilibrium Bondi solution for any 
magnitude of radial magnetic field. The magnetized Bondi problem is set up in the Schwarzschild 
metric (a = 0). The radial grid extends from just outside the event horizon at r = 2.20 M to 
r = 25.0 M; the critical point is at Vcrit = 8.0. Grids of 64, 128, 256, 512, and 1024 points are used, 
and the magnetic field at the critical point is set to /3 = 100, 10, and 1, where the /? parameter 
sets the grid-averaged ratio of gas to magnetic pressure. The Bondi solution is generated in the 
initialization routine by a numerical root-finder on a logarithmically scaled grid of 1024 points. 
This solution is then sampled at the required resolution for the particular test run (to ensure 
proper alignment of the coarser grids). To complete the test, the unmagnetic case is also presented. 
The solution is allowed to evolve to a time sufficient to allow numerical convergence on each grid 

(t final = 100 M). 

Figure 10 shows the results of spatial accuracy tests. The spatial accuracy of variable X is 
obtained from the Li norm over an interval rmin < ^ < fmax, 

Li{X) = [ WXfinai - XoWdr, (62) 

where Xfinai denotes the converged solution, and the original solution, as described above. 
Different regions of the flow arc dominated by numerical errors cither from a particular routine in 
the numerical solver or a particular code variable. Here, the bounds of integration were chosen to 
range from r^m = 7.0 to rmax = 25.0, which represents a region of the flow where pressure error 
dominates. Figure 10 is for the variables d, and is typical of error curves for all code variables. The 
code is clearly first-order accurate for these variables. 

In the split monopole test, we initialize a 2D (r, 6) grid with the radial Bondi solution, but 
the radial magnetic field has opposite polarity above and below the equator, hence establishing a 
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Fig. 10. — Plot of spatial error for Bondi test, for values of /? = 100 (crosses), 10 (diamonds), 1 
(triangles) and the unmagnetized case (squares) for variable d . 

current sheet along the equator. This test is made into a numerical stability test by applying small 
enthalpy perturbations to the initial state (without this, the code is simply a set of ID radial tests). 
The test is repeated with increasingly large values of magnetic field, until the code breaks. There 
is no known analytic solution for the evolution of a current sheet, but it can be used as a probe of 
the code's stability in the presence of strong, fluctuating current sheets. The split monopole test 
was carried out with (3 = 100, 10, 1 at the critical point, and a grid of 180^ scaled logarithmically 
in both the r and 9 coordinates, with an initial 1% random enthalpy perturbation. The test was 
allowed to run to t = 1250M (10 free-fall times). For /3 = 100, 10, the code variables show no 
visible change over the course of the run; the enthalpy perturbations do not grow, nor do they 
trigger a numerical instability. For (3 = 1, the late-time solution shows evidence of reconnection at 
the equator near the horizon, and erratic behavior. The effective /3 where reconnection occurs is 
of order 10^^. Given the moderate resolution, the results are consistent with the spatial accuracy 
tests in that the region near the horizon requires high resolution in order to minimize numerical 
error. 
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5.2. Magnetized Gammie Inflow 

Gammie (1999) describes a simplified ID model for inflow from the inner edge of an accretion 
disk in the Kerr geometry that serves as an excellent test problem for relativistic MHD. This 
solution is similar to the Weber & Davis (1967) rotating magnetized wind which was used as a 
test problem Stone et al. (1992). The model assumes that the inflow is cold {h k, 1), time- 
steady, purely radial, confined to the equatorial plane, and initiated from the edge of an accretion 
disk located at the marginally stable orbit {rrnso 

). At 

Tmso the solution is assumed to satisfy the 

conditions V {vmso 

) = and Of = V'i>{rmso) = l/(r mso + ci) • The flow is characterized by a fast 
magnetosonic critical point. The numerical procedure to find a particular solution to this inflow is 
described in Appendix D. 

Our test consists in replicating the solution described by Gammie et al. (2002) for a flow 
near a Kerr black hole with a = 0.5 {r horizon = 1-866, rmso = 4.233, Clp = 0.108588), with 
= -0.5, Fm = -1, Fl = -2.815344, Fe = -0.9083782. The critical point for this test case is 
{rcriuUl^it) = (3.616655,-0.04054696). 

Figure 11 shows a plot of the solution for this test problem on a logarithmically scaled grid 
of 1024 points. The inner edge of the grid is at Vmin = 2.00 and the outer edge of the grid is at 
Tmax = 4.04. In Boyer-Lindquist coordinates, the event horizon cannot be part of the computational 
domain because of the coordinate singularity. In the Gammie inflow test, the marginally stable 
orbit, which is the formal outer boundary of the problem, is also excluded since the density is 
divergent as r approaches Vmso from below. The plot overlays the late-time solution (after 1000 
time steps) on the initial state, but the two curves are indistinguishable at the chosen plot scale. 
Figure 12 shows the results of convergence tests carried out analogously to those for the Bondi 
problem, with the range of integration going from Vmin = 2.7 to r^ax = 4.04. Here we have an 
example of a highly-resolved near-horizon region, and the convergence rate is close to second order. 
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Fig. 11. — Plot of four- velocity components, density, and magnetic pressure for the Magnetized 
Gammie Inflow test, with overlaid plots for i = and t = 15M. Grid 1024 points. 
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Fig. 12. — Plot of spatial error for Magnetized Gammie Inflow, (left) variable d and (right) variable 
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5.3. Constant-Z Disks 

The final set of tests involve constant-/ disks (/ = —U,j)/Ut, the specific angular momentum), 
solutions of the axisymmetric GR hydrodynamic equations described in HSWa. Here, we add 
weak poloidal magnetic field loops that overlay the hydrodynamic solution; this will trigger the 
MRI. We quantify the strength of the magnetic field through the /3-parameter, the ratio of the 
volume-averaged gas pressure to magnetic pressure. 

The initial magnetic field is obtained from the definition of F^i, in terms of the 4-vector 
potential, A^, 

F^u = d^A, - d^A^. (63) 
Restricting the field to poloidal loops is done by using A^ = {At, 0, 0, A^), where 

I for p < pmin 

Using (63), it follows that B'^ = —dgArj) and = 5^^^. This choice of A^p produces poloidal field 
loops that coincide with isodensity contours. The constant k is set by the input parameter /3. The 
constant pmin sets a suitable minimum density within the disk, and it is chosen to keep the initial 
magnetic field away from the outer edge of the disk. 

An axisymmetric (2D) magnetized constant-Z torus was set up in both the Schwarzschild and 
Kerr metrics. For the Schwarzschild case, the disk has a specific angular momentum I = 4.5, 
an initial pressure maximum at r = 15. 3M, and an orbital period at the pressure maximum of 
Torb = 376M. For the Kerr case, the black hole is maximally rotating, a = 1, and the disk has 
a prograde specific angular momentum I = 4.3, an initial pressure maximum at r = 15.4M, and 
an orbital period at the pressure maximum of Torb = 386M. This particular choice of parameters 
yields a disk that is similar to the Schwarzschild case, although the equipotentials which define the 
overall disk structure are notably different, as can be seen in the left column of Figure 13. For both 
tests the average field strength is /? = 100, the initial state is perturbed with random 1% enthalpy 
fluctuations. The simulations are run for 10 orbits. The MRI develops after a few orbits (center 
column. Fig. 13), and is soon fully developed (right column. Fig. 13). By the tenth orbit, the 
MRI-induced turbulence has settled down considerably, as is to be expected since the imposition 
of axisymmetry precludes the development of the azimuthal modes which sustain the MRI. The 
observed qualitative differences between the Schwarzschild and Kerr cases can be attributed to the 
different shapes of the potential cusp near the horizon. 

The density plots bear a strong resemblance to the previous axisymmetric MRI studies of 
Hawley (2000) in a pseudo-Newtonian potential; this indicates that the numerical solver is well 
behaved in the moderately strong gravitational fields that exist near the pressure maximum, where 
the implied code comparison is taking place. Since MRI studies are the main area of application 
of this code, it is reassuring that we can both trigger the MRI in a general relativistic treatment 
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of MHD, and that we qualitatively reproduce established results. Gammie et al. (2002) obtain 
similar results with their MRI torus test. 




10 20 30 40 10 20 30 40 10 20 30 40 




10 20 30 40 10 20 30 40 10 20 30 40 



Fig. 13. — Plot of log density for the Schwarzschild (top row) and Kerr (bottom row) constant-/ 
disks. The left column is the initial state, the middle column is from t = 1.5 orbits, and the right 
column corresponds to t = 3.0 orbits at the initial pressure maximum. 
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6. Discussion 

In this paper wc document the development of a General Relativistic MHD code that is a direct 
extension of the general relativistic hydrodynamic solver of Hawley, Smarr and Wilson (1984b). As 
before, the equations of relativistic MHD are written in a nonconservative form, supplemented here 
by the relativistic induction cqiiation. The magnetic field is added to the definition of the inertia 
for the four momentum and the magnetic field adds several source terms to its evolution. 

The equations axe evolved on a staggered grid using the time-explicit model of "transport + 

source." The induction equation is evolved using the Constrained Transport structure of EH88. We 
implement a version of CT that is inspired by the Method of Characteristics technique employed in 
the popular ZEUS code (Stone &: Norman 1992), yet is simpler to implement in a general relativistic 
framework. 

We have performed a number of ID and 2D tests to validate the code and to note its limitations. 
The tests include a suite of non-relativistic and relativistic MHD shock tube tests, the fast and 
slow rarefaction and slow magnetosonic shock tests of Komissarov (1999), and simple Alfven pulse 
propagation. The formulation for these test and their results are documented to establish a basis 
for comparison against other existing and future codes. 

The main shortcoming of our code revealed by these tests is due to the already-known limita- 
tions of the artificial viscosity algorithm of the underlying hydrodynamic solver. Improvements to 
this algorithm are an ongoing area of research, but is important to stress that these limitations do 
not pose an undue restriction on our main area of study: dynamical, MHD-driven accretion flows 
around spinning black holes. 

Although the set of general relativistic black hole accretion test problems is limited, we have 
also verified that the code reproduces the standard results for Bondi flow in the presence of magnetic 
fields, and have observed the triggering and evolution of the MRI in thick accretion tori for both 
the Kerr and Schwarzschild black holes. The MRI evolution can be compared qualitatively with 
pseudo-Newtonian simulations (e.g. Hawley 2000). One new test problem with an analytic solution 
is the magnetized Gammie (1999) inflow test. This solution examines the behavior of a dynamically 
important magnetic field in the strong-field region of the Kerr metric. The results of these tests 
provide evidence that this new solver is correctly reproducing magnetized flows outside the event 
horizon of a Kerr black hole. This new GR MHD code should provide a valuable tool with which 
to study accretion flows in the Kerr metric. 

This work was supported by NSF grant AST-0070979 and NASA grant NAG5-9266. We thank 
Charles Gammie with whom we collaborated to develop and apply a suite of code tests for general 
relativistic MHD. 
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A. Proof or Various Results 

Some of the following proofs use two identities for contractions of the Levi-Civita tensor, 

a/357 , — _xl3S'r ( 

e'''''^e^,3^. = -2S'\,, (A2) 

and identities for the 5-symbol, 

S'^^, = S^S^^-SX, (A3) 

To prove (10), it is easiest to work backwards from the result. Expand the right-hand side of 
(10) using (8) and (6); the resulting expression simplifies through identities (Al) and (A4), followed 
by the velocity normalization condition, infinite conductivity, and antisymmetry of F^'^: 

= —(UpU^Fi,,- Fp, + F,fs - F^p + u,U'^Fp,- Up F,^) 

= F 

To prove (12), expand the EM component of tensor (7) using (10), and simphfy using identities 
(A1)-(A4), followed by the velocity normalization condition, and orthogonality condition (11): 



■'-[EM) 



1 

47r 

1 



g'^^S's^^B'U'BpUa 
1 / gl'-B'^Bo, 



J s:Pcr T30 Tre 



4^>. ^5e,^ - ^,u.-—6^s:b'U^B,U. 



^Bi^B" - 5" Ba V U" 



To prove (15), the second of Maxwell's equations, Vf/F'^" = 0, can be rewritten by substitut- 
ing definitions, and using identities (A2) and (A3): 

= -VaS'^^bPU'', 
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which simplifies, using (A3), to Va (C/" - 6" U^) = 0. 

To prove (38), we begin with an identity derived from the expression for W (18), 



a-2 = = -{gtt + 2gt4, + grr {V'f + gee {V'f + g^^ {V^) 



2.-1 



or 



Now expand 



and substitute 



gtt + grr {vn^ + gee {v'T + g^ (^'^) = - '^9t4> V^. 



gttib'f + grr {b'f + gee (b'f + 5<a<a (b'^f + '^gtct> b' b\ 



(A5) 
(A6) 

(A7) 
(A8) 



for the spatial components of the magnetic field. Making use of the above identity, it follows that 



{b'f + 2gt^ 



4itW ^ 



+ b*V'^ 



(A9) 



+ 



+ 



AttW ^ 
1 



gr 



B'V + geeB<^V'^ +g^^B't>V't> 



grr{Bn^+gee{Br + g4,4,{B^) 



Now substitute the expression for 6* (37) and simplify, noting that 

.V'f' + gtd> = 



gtt — gt't> 

g4><pv- + gt4> = . ,.,2 



gtt g(j)(j) _ {^gi't>y 



(AlO) 



to get, after a few steps, expression (38). 



To derive (51) we first obtain the components of the fluid-frame magnetic field, which are, to 
leading order. 



6* 
b^ 



g. 



-y/47r7' 
WB"" 



By w v yy b^ 

i~ gxx 



B^ 



gx 



and also note the approximate value for W, 



(All) 
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We then expand the momentum equation and substitute the expressions for by and 6*. The left-hand 
side of the momentum equation reads 



dt {Sy -aby 6*) = dt {gyy {ph+ ||6f ) VF^ - gyy a by b' 



dt [gyyph — V^+gyyWbf — Vy 



9xx {V' f 



gxx9yyB-Byv- 



gyy {phW^ + \\b^ _ g..gyyB-V- ^^^^ _ 



" a (•\/47r7)' 

Recall that x stands for either r or 9, and y stands for a transverse direction, either r, 6, or ^ 

As for the left hand side, we neglect the purely geometric terms and retain only 

-^d,^{SyV--abyb-) = -±d,^ (gyyiph+\\bf)^VyV--agyybyb'''^ 

^ dxV^ \ 9yyPh — VyV -agyy^^^—^] 



2 ^/x 



_gyyphW V-^^ ^ agyyB- 



Combine these results and replace dt {By) with the results for the induction equation, 

{phw^+\\br) 



9xxB^V^ 


A'K^{phW^ + II 








47r7(/9/iVF2 + II 


^11') 


gxxB^ 




4.T:-i{phW^ + \\ 




a^B"^ 




ATr^{phW'^ + \\ 


5'f) 



+ z — . , . u,u2. 9x {By) 

{vy)-v'' d,{By)) 

+ Z . . u,u2. 9x{By). 

By regrouping these terms and substituting the approximate expressions for W and ||6|| given 
above, equation (51) follows immediately. A similar procedure yields expression (56) for transverse 
velocities when x = (f). 

To prove that the Bondi solution is unchanged by the presence of a radial magnetic field, 
we revisit the calculations presented in HSWa. We are looking for a time-independent purely 
radial solution to the equations of GRMHD in the presence of a constant CT radial magnetic field 
B'^' = Fg^ = B in the Schwarzschild metric. We will work with the more general formalism of 

equations (2) and (3). 

The equation of continuity (2) is not dependent on magnetic fields, so 

V^(pC/'^) = 0^V,(/9C/^) = 0, (A12) 
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leads to the same conserved quantity as in the hydrodynamic case, 

^/^gpU'^ = Ci. (A13) 

The energy-flux conservation equation is the t-component of (3), 

V^T'^t = ^ VrT't = 0, (A14) 

and contains an expHcit dependence on the magnetic field. However, by expanding the r — t 
component of the energy momentum tensor, 

T't={ph+\\bf) U'Ut-b'bt, (A15) 

and using equations (32) and (37) to obtain the fluid-frame magnetic field. 



we obtain after simple substitutions 

rt={ph+\\bf) U'Ut-\\bfU'''Ut. (A17) 

Clearly, the magnetic contributions cancel (this is simply restating the general result that there is 
no net force due to magnetic fields parallel to the direction of motion). We therefore recover the 
hydrodynamic expression, 

VrT''t = {ph)U'-Ut = 0, (A18) 
which leads to the same conserved quantities as in the hydrodynamic case, 

^pU'-hUt = C2, (A19) 

and 

hUt = C3. (A20) 



B. Amplitude of Alfven Pulses 



To prove (61), we used the method of characteristics. We write the momentum and induction 
equations as a matrix equation, 

Ad^U + BdtU = (Bl) 
where U{x,t) = {Vy{x,t),By{x,t)f and 

«2 



A 



_jgx yx 



B 



1 X 
1 



B"" (1 + rf 



(B2) 
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subject to the initial state Uo{x,0) = (/(x),0)'^ We obtain the eigenvalues from the roots of 
F{X) = detA — XB, which are simply the expressions for va^'^\ (59). 

To solve this system, we seek a matrix P such that PBU = z where Zi{x,t) = gi{x — Xit) 
and g{x, 0) = P B Uo{x,0). The rows of matrix P satisfy — Aj B^ Pi = 0. It is straightforward 
to show that the matrix P has the form, 



We obtain the results 



c 

L-C 



(B3) 



which yield 



Cy^ + (i + Cx)H^ = Cfix-VA't), 
-Cy^ + (i-Cx)H^ = -c,!{x-VA^t\ 



1-Cx 



j(x -VA t) + 



1 + Cx 



f{x-VA'^t), 



\ {f{x-VA-t)-fix-VA^t)) 



(B4) 
(B5) 



(B6) 
(B7) 



As discussed in the text, the velocity amplitudes are clearly asymmetrical, except for the case where 
= 0, and the magnetic field amplitudes are always equal, but opposite in sign. It is easy to 
show that the ratio of the velocity amplitudes is indeed the ratio of the gamma factors, W{va^^^)- 
It is also possible to verify by direct computation that the amplitudes of the pulses given in the 
text agree with the above formulas. 



C. MHD Shocks 

Relativistic MHD shocks have been discussed by Anile (1989) and Komissarov (1999). Their 
work builds upon the monograph of Lichnerowicz (1967) who systematically developed the rela- 
tions for GR MHD shock invariants. We summarize here the relevant results from Lichnerowicz, 
translating symbols to our own notation. We use these results to construct our shock tests. 

MHD shocks are discontinuities in one or more of the MHD variables whose propagation 
through spacetime lies on a hypersurface S. It is useful in the analysis of these shocks to construct 
a 4- vector normal to S. Shocks can be characterized by invariant scalar and 4- vector quantities — 
the jump conditions — which are derived readily from the equations of GRMHD. MHD shocks can 
be assigned to one of two categories, tangential and non-tangential, and this classification hinges on 
whether the invariants have non-zero projections onto the normal n^. Slow and fast magnetosonic 
shocks are examples of non-tangential shocks, so it is this type of shock that is of interest here. 

The analysis of non-tangential shocks was shown by Lichnerowicz to reduce to the study of 
five scalar quantities, p, h, rj = h^n^, A = pW^Un, and which are related through five 



positive-definite scalar invariants (derivable from the jump conditions), 



[A] = 0, 

[B] = [hv]=0, 



(CI) 
(C2) 




= 0, 



(C3) 



(C4) 



(C5) 



where [X] = — X- indicates the change in a quantity X across the shock, and 




}?\\-A^H, 



(C6) 
(C7) 



(C8) 



Magnetosonic shocks arc characterized by the parameter a. The determination as to whether a 
magnetosonic wave is fast or slow hinges on the sign of the pre- and post-shock values of parameter a. 
Physically realizable shocks are characterized by a(+) a(_) > 0. A slow shock has q;(_|_) < a(_) < 



To construct an initial state, we solve the scalar invariants using the Mathematica FindRoot 
function. First, we postulate a left state in the rest frame of the shock, where the normal vector 

n'^ has an especially simple expression (n'^ = (0, 1,0,0)). Next, we set up the scalar invariants as a 
system to be solved by the FindRoot function, supplying an initial guess for the right state. This 
is an iterative procedure, requiring adjustments to the initial guess. Once a solution is found, it is 
possible to construct a moving shock by repeating the above process with a small non-zero shock 
speed, for which the normal vector now has the form = W{vsh) {vsh, 1; 0, 0). Once a solution is 
found the process is repeated, gradually increasing the shock speed and using the previous converged 
right state as a starting point for the solution at the new shock speed. This is a pragmatic solution 
to a problem that has a high degree of complexity, and it allows us to set up a family of closely 
related shock solutions to aid in testing. 

Given that the FindRoot function has converged on a right state, it is possible to compute the 
parameter a, and classify the shock. It is also prudent to verify that the left and right states are 
true solutions by verifying that they satisfy the shock invariants by direct computation. Once this 
is done, the left and right states can be transferred directly to a numerical grid. 



and a fast shock has < q;(_|.) < «(_). 
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D. Determining Conserved Quantities for Magnetized Gammie Inflow 

The magnetized Gammie inflow is characterized by three conserved quantities, 



rr m 



Fe = -2Trr'^T\ = -2Trr 



rr m 

pU'Ut + ^-f-FroFte 



(Dl) 
(D2) 

(D3) 



which are readily derived from the equations of baryon and momentum conservation. Note that 
Gammie writes V = r'^ g^^ g^^. From the induction equation, 



drB'' = 0, 
VB't'^ = 0, 



(D4) 



we obtain = —Fq^ = const, and Fiq = B^ — V'^B'^ = Qp Fq^ which follows from the boundary 
condition. Infinite conductivity {V^ F^y = 0) yields 



Fre='-^{u'»-^pU'). 
Finally, the velocity normalization condition yields 

-1 = 5** iUtf + 2g"t' Ut + /"^ iU^f + iU'f/g'' 
where the mix of upper and lower indices was made to coincide with code variables. 



(D5) 



(D6) 



The critical point is a saddle point in the two-dimensional parameter space of the energy 
function Fe = FE{r,U'^;Fi,). It is apparent that Fe as given above is not a simple function of 
these parameters, and it is necessary to substitute the other expressions in order to obtain the 
desired form for Fe- Since the resulting expression for Fe is quite lengthy, we will only sketch 
the procedure. The process of substitution and elimination described here was carried out using 
Mathematica. 

First, the conserved quantity Fm is used to obtain p{r,U^) = Fm / {2 tt U^) . Next, the 
conserved quantity Fl is used to isolate Us, 



Us 



Fl 



+ 



V 



2Trr'^U'' 47rr2 



Fed) Fr0 



1 



pU^ 



(D7) 



Now substitute the expressions for p{r, W) and FrO expressed in terms of U^ and Ut, to yield an 
expression of the form U^ = A + B Ut where A and B are expressions involving the constants Fl, 
Fm, Fg^j,, and VI p, as well as U^ and metric terms. Now substitute this expression for U^ in the 
velocity normalization condition and "solve" the resulting quadratic for Ut, denote this solution 
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(Ut) to emphasize that there are two distinct roots. Finally, take the conserved quantity Fe 
and substitute the expression for F^g (expressed in terms of and Ut) and then perform a final 
substitution for and Ut by (Ut)^. Simplify the resulting expression, which has the desired 
functional form: (r, [/'' ; F/^ ) . 

The critical point must satisfy the conditions 

drFE{rcrit,U^rit'^FL) = 0, 

durFE{rcrit,Ul^it\FL) = 0, 
and the conserved quantity Fe must also satisfy the constraint at the outer boundary 

FE{rcrit,U^rivFL) = FE{rmso,0; Fl)- 

These three conditions allow us to solve for rcrit, t^crif Fl- 

In order to find the critical point, wc will eventually use the FindRoot function of Mathematica, 
but we first need to obtain a rough location for the saddle point. To do this, we produce a coarse 
contour plot over the range Vfiorizon < r < Vmso and a large range of U^, say —0.5 < ?7'' < for a 
given choice of root of {Ut)"^. If the contour plot does not show evidence of a saddle point, we choose 
the other root. Figure 14 illustrates this for our specific test problem. The saddle point is shown 
to lie near (vcriti ^Irii) ^ (3.65, —0.04). This is then used as an initial guess for the Mathematica 
FindRoot function, from which the parameters used in the test were obtained. 




3.2 3.4 3.6 3.8 4 4.2 



Fig. 14. — Contour plot of the function FE{r, U'^;Fi,). The saddle point is located at the intersection 
of the dashed contours. The ingoing solution follows the steeper of the dashed contours; the outgoing 
solution follows the other dashed contour. 

To generate an initial state within our code, we apply a numerical root-finder to the energy 
function FE{r,U^;FL) (coded using the Mathematica FortranForm command), and generate a 
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numerical solution for (using the parameters in the text). Once is found, we construct Ut 

and U(f, using the algebraic expressions obtained from Mathcmatica (also coded using FortranForm) . 
Once these variables have been initialized the remainder of the code variables can be obtained in a 
straightforward manner. 
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